;;------ solve non-linear eqns by newton-sceant method

;; Problem f(x) = 0, x is in R^n

;; load Least Squares method to solve z_k
(load "/home/cl/gits/scheme/numerical-analysis/linear-eqns/least-squares.scm")
;; stream pkg
(load "/home/cl/gits/libs-for-chezscheme/stream.ss")
(import (stream))
;; matrix-transport: a matrix is a vector of vectors(rows)
(define (matrix-transport A)
  (let ([n (vector-length A)])
    (vector-map (lambda (i) (vector-map (lambda (row) (vector-ref row i))
                                        A))
                (vector-range n))))

;;next approximate solution
(define (nextsol x0 x1 f n)
  (let* ([H (vector-map - x0 x1)];H is a diag-matrix
         [xi-func (lambda (j); xi-func(j) => xi_j
                    (vector-map (lambda (l)
                                  (if (= l j)
                                      (+ (vector-ref x1 l)
                                         (vector-ref H j))
                                      (vector-ref x1 l)))
                                (vector-range n)))]
         [Gamma (matrix-transport
                 (vector-map (lambda (j)
                               (vector-map - (f (xi-func j)) (f x1)))
                             (vector-range n)))]
         [z (Least-Squares-method Gamma (f x1))])
    (vector-map - x1 (vector-map * H z))))

;; newton-sceant method
(define (newton-sceant-stream f x0 x1 n);;need 2 initial points
  (stream-cons x0
               (newton-sceant-stream f x1 (nextsol x0 x1 f n) n)))

#| ;;test
(let* ([f (lambda (x)
            (vector (- (* 3 (vector-ref x 0))
                       (cos (* (vector-ref x 1) (vector-ref x 2)))
                       0.5)
                    (+ (expt (vector-ref x 0) 2)
                       (* -81 (expt (+ 0.1 (vector-ref x 1)) 2))
                       (sin (vector-ref x 2))
                       1.06)
                    (+ (exp (* -1 (vector-ref x 0) (vector-ref x 1)))
                       (* 20 (vector-ref x 2))
                       (* 10/3 3.1415926)
                       -1)))]
       [solutions (newton-sceant-stream f
                                        (vector 0. 0. 0.)
                                        (vector 0.01 0.01 0.01)
                                        3)]
       [steps 8])
  (display "------ f(x) -------\n")
  (for-each (lambda (i)
              (display (f (stream-ref solutions i)))
              (newline))
            (range steps))
  (display "-------- x --------\n")
  (for-each (lambda (i)
              (display (stream-ref solutions i))
              (newline))
            (range steps))
)
|#






